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Abstract 

Molecular modelling and simulation as well as four equations of state (EOS) are ap- 
plied to natural gas mixtures regarding Joule-Thomson (JT) inversion. JT inversion 
curves are determined by molecular simulation for six different natural gas mixtures 
consisting of methane, nitrogen, carbon dioxide and ethane. These components are 
also regarded as pure fluids, leading to a total of ten studied systems. The results 
are compared to four advanced mixture EOS: DDMIX, SUPERTRAPP, BACKONE 
and the recent GERG-2004 Wide- Range Reference EOS. It is found that molecu- 
lar simulation is competitive with state-of-the-art EOS in predicting JT inversion 
curves. The molecular based approaches (simulation and BACKONE) are superior 
to DDMIX and SUERTRAPP. 
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1 Introduction 



Due to the eminent importance of natural gas, knowledge on its thermodynamic 
behaviour and appropriate property models are of great interest. A property often 
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needed for applications of natural gases is the adiabatic (or isenthalpic) Joule- 
Thomson ( JT) coefficient which is defined as the derivative of the temperature 
T with respect to the pressure p at constant enthalpy h and constant composition 

X 



where Cp is the isobaric heat capacity. The JT inversion curve, connecting all state 
points with n—O, divides the pressure-temperature plane into two regions. In the 
lower region, /x is positive so that an adiabatic expansion leads to a decrease in 
temperature. In the upper region, /i is negative. It can be shown that the cooling 
effect is maximized if an expansion starts from the inversion pressure. 

The experimental determination of a fluid's JT inversion cocflicicnt demands 
precise measurements of volumetric and caloric properties, while the JT inver- 
sion curve extends over a broad range of temperature and pressure. Temperature 
and pressure reach up to a five-fold and twelve-fold of their critical values, re- 
spectively. Experimental JT inversion curve data are therefore scarce, sometimes 
unreliable [1], and mostly available only for pure fluids [2]. An example for mix- 
ture data is the work of Charnley et al. [3] on e.g. carbon dioxide -|- nitrous oxide, 
carbon dioxide -|- ethylene. A good overview over experimental data is given in 




(1) 



or, using basic thermodynamics relations. 





[2]. 
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For industrial applications there is a need for a proper representation and also for 
the prediction of JT data, as it is not feasible to measure it for all relevant and 
often very different blends. Some authors, e.g.. Miller [4] or Gunn et al. [5], have 
proposed direct representations of the JT inversion curve, which are correlations 
in terms of reduced temperature and pressure. They provide rule-of-thumb data 
for simple fluids, but have little predictive power. Signiflcantly more valuable 
are proper equations of state (EOS) that contain much more thermodynamic 
information and are valid for a broad range of state points. Thus, extensive efforts 
are made to use EOS for predictions of the JT inversion curve. Examples for the 
use of cubic EOS are modified versions of Peng-Robinson [6,7], Redhch-Kwong 
[6,8], Soave-Redlich-Kwong [7], Patel-Teja [6] or other cubic EOS [9] with varying 
parameter functions. Both type of cubic EOS and type of parameter function 
strongly influence the JT inversion curve, particularly in the high temperature 
region. A given combination might yield good results for a speciflc fluid, but fails 
for others [6,7,8,9]. Hence, it can be concluded that cubic EOS are not generally 
reliable for JT inversion curve predictions. 

Among mixtures, natural gases are the ones that were investigated most exten- 
sively both experimentally and theoretically so that very reliable thermodynamic 
data and models are available. Therefore, natural gases are excellent test systems 
to validate thermodynamic models for mixtures. 

Natural gas from the rig is a mixture of typically seventeen components (con- 
taining methane, nitrogen, carbon dioxide, ethane, propane) [10], but usually its 
main component is methane [11,12]. As a natural product, it has a great vari- 
ability in composition and, depending on conditions in the formation process, 
considerable quantities of nitrogen (up to 60 mole %), carbon dioxide (up to 
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mole 50 %) or ethane (up to mole 20 %) are encountered [11]. 

For a number of pure natural gas components, reference EOS have been devel- 
oped based on a vast experimental data set considering different thermodynamic 
properties, e.g., methane [13], nitrogen [14], carbon dioxide [15] and ethane [16]. 
Reference EOS have an empirical background, but they are parameterized ex- 
tremely carefully, taking also available experimental JT coefficients into account. 
Hence, they are regarded in this work as the best available information. 

For mixtures, the National Institute of Standards and Technology (NIST) [17] 
provided two classical phenomenological EOS, i.e. DDMIX [18,19] and SUPER- 
TRAPP [20]. DDMIX is an implementation of the NIST extended corresponding 
states model for mixtures, whereas SUPERTRAPP is based on both a modi- 
fied Peng-Robinson EOS and the NIST extended corresponding states model for 
mixtures. Both were parameterized to experimental pure substance and mixture 
data. Particularly SUPERTRAPP is often used in the literature as a property 
model for designing cooling cycles with mixed coolants, e.g. [21,22]. 

There are also physically based EOS that take the different molecular inter- 
actions, like dispersion or polarity, explicitly into account; an example is the 
BACKONE-EOS [23]. Such EOS can be parameterized for real substances with 
a very small experimental data set, e.g. a few vapour-liquid equilibrium data 
points, as they have a good predictive power. Furthermore, the GERG-2004 
Wide-Range Reference EOS [24] has become available recently, which applies 
the concept of reference EOS to mixtures. Here too, a vast set of experimental 
data, both pure substance and mixture, was used for the development. 
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Molecular modelling and simulation offers an interesting alternative approach for 
predicting thermodynamic properties. Instead of describing those macroscopic 
properties directly, the intermolecular interactions are described. Previous work 
from our group [25] has demonstrated the good predictive power of molecular 
models for JT inversion curves for different pure fluids, but also for the mixture 
air which has been modelled as a ternary system containing nitrogen, oxygen and 
argon. The objective of the present paper is to validate the predictive power of 
available molecular models by comparing the results to the four advanced EOS 
mentioned above using different natural gas mixtures and their most important 
pure components regarding JT inversion. 

2 Molecular model and simulation method 

Most publications on JT inversion curves by molecular simulation are based on 
the spherical Lennard- Jones (LJ) model [26,27,28,29,30,31,32,33], which is ap- 
propriate only for simple molecules like methane and the noble gases. However, 
this simple potential model is well suited to further develop molecular simulation 
techniques for determining JT inversion as reliable simulation data is available 
for comparison. E.g., Colina et al. [34] have chosen two different routes, i.e. 
via compressibility [30,31] or via thermal expansivity [32], to simulate the JT 
inversion curve. Work on more complex fluids is still scarce, examples are Es- 
cobedo et al. [33] (nitrogen), Chacfn et al. [35] (carbon dioxide), Lisal et al. [36] 
(R32), Kristof et al. [37] (hydrogen sulphide), or a recent work of our group [25] 
dealing with fifteen different pure substances (including methane, carbon diox- 
ide, R134a, R143a, R152a) and the mixture air. Two publications, dealing with 
multi-component natural gas mixtures, should be mentioned: Escobedo et al. [33] 
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predicted the JT inversion curve for a seven-component system, but compared it 
to cubic EOS data only. A favourable comparison between simulation and exper- 
iment was presented by Lagache et al. [38] regarding the JT inversion pressure 
of a 20-component mixture for one specified temperature. 

The effective 2CLJQ pair potential was used as molecular model here to describe 
the intermolecular interactions in all cases. This can be done, as only the four 
most common components methane, nitrogen, carbon dioxide and ethane were 
considered, being well suited for this modelling approach. The 2CLJQ potential 
is composed out of two identical Lennard- Jones sites a distance L apart (2CLJ) 
and a point quadrupole with momentum Q placed in the geometric centre of the 
molecule oriented along the molecular axis 

^i2CLjQ(r-ij, Wj, ujj, L, Q) = ii2CLj(nj, Wj, u?j, L) + uq{rij, Wj, w^, Q), (3) 

wherein M2clj is the contribution of the four Lennard- Jones interactions 

2 2 
a=l 6=1 

The quadrupolar contribution uq, is given by [39] 

3 

Herein, r^j is the ccntrc-ccntrc distance vector of two molecules i and j. Tab is one 
of the four Lennard- Jones site-site distances; a counts the two sites of molecule 
i, h counts those of molecule j. The vectors and ujj represent the orientations 
of the two molecules i and j. /q is a trigonometric function depending on these 
molecular orientations, cf . [39] . The Lennard- Jones parameters a and e represent 
size and energy, respectively. 2CLJQ models have the four state independent 



-) 
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a 

Tab 



(4) 
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model parameters o", e, L, and Q, which have been adjusted to experimental 
vapour pressure, bubble density and critical temperature in a recent work of our 
group [40]. All pure substance parameters are given in Table 1. The spherical 
non-polar LJ model for methane is a limiting case of 2CLJQ models, where L — 
and g = 0. 

To perform simulations of mixtures, a molecular mixture model is needed. On 
the basis of pairwise additive pure fluid models, molecular modelling of mixtures 
reduces to modelling interactions between unlike molecules. Here, the modified 
Lorentz-Berthelot combining rule with one adjustable binary interaction param- 
eter ^ was used for each unlike Lennard- Jones interaction 



The state independent parameter ^ was adjusted to one experimental vapour 
pressure of each binary mixture in prior work of our group [41,42,43]. Table 2 
reports the six binary parameters of the quaternary natural gas mixture model. 
The unlike quadrupolar interactions are treated in a physically straightforward 
way, following the laws of electrostatics without any binary parameters. It should 
be pointed out, that exclusively experimental VLE data were used in the param- 
eterization of the molecular mixture model, but no caloric data. 

For the calculation of JT inversion curves on the basis of a given model by molecu- 
lar simulation, several methods have been proposed in the literature [30,31,32,44]. 
Here, as in prior work [25], an intuitive and straightforward method was used. 
To calculate one JT inversion pressure p for a given temperature T, a series of 



CTj + tTj 

2 



(6) 



(7) 
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simulations, generally from 5 to 10, were made around the expected result, cov- 
ering typically a rather large pressure range of 20 MPa. In these simulations, 
the enthalpy h and its partial derivative with respect to pressure {dh/dp)T,x 
were calculated. Both data sets were fitted simultaneously by a second order 
polynomial vs. pressure at that particular temperature. The inversion pressure 
corresponds simply to the minimum value of enthalpy, i.e. the minimum of that 
quadratic fit. 

Molecular dynamics simulations were performed in the isobaric-isothermal (NpT) 
ensemble, using Andersen's barostat [45] and isokinetic velocity scaling [46] for 
thermostating. After 6 000 equilibration time steps, the residual enthalpy [47] 

I ( N N ] 

h'"' = T7- < < E E U2CLjQ{rij,uJi, ujj, L,Q)> +p<V > \ - /cbT, (8) 
-'^ [ 1=1 j>i J 

was averaged over 200 000 time steps, where the first term indicates the simula- 
tion average over the intermolecular potential energy and < y > is the average 
of the extensive volume. The partial derivative was obtained by a fluctuation 
expression [48] 

where H^'^^ is the extensive residual enthalpy. The ideal part of the enthalpy is 
not relevant, as it is not pressure dependent. 

Depending on the density of the state point, the membrane mass parameter M 
of Andersen's barostat [45] was chosen from 10~^° to 10"^^ kg/m^. The inter- 
molecular interactions were evaluated explicitly up to a cut-off radius of Tc = 5a 
and standard long range corrections were used for the Lennard- Jones interaction. 
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employing angle averaging as proposed by Lustig [49] . Long-range corrections for 
the quadrupolar interaction are not needed since they disappear. A total number 
of = 1372 molecules were initially placed in a fee lattice configuration into a 
cubic simulation box, where the density of the system was chosen close to that 
expected from an EOS, if available, otherwise estimates were used. 

3 Results 

JT inversion curves are compared for six different systematically chosen gas mix- 
tures consisting of the four main natural gas components methane, nitrogen, 
carbon dioxide and ethane. Firstly, the three methane containing binary sys- 
tems that can be formed from these four components were studied at equimolar 
composition, where the highest effect of mixing can be expected. Following this, 
secondly, all three ternary systems containing methane were studied, again at 
equimolar composition, i.e. with a mole fraction Xi — 1/3 throughout. Finally, 
also the four pure fluids were investigated. 

Simulation results for pure methane and carbon dioxide were taken from prior 
work [25], but for nitrogen and ethane as well as for all mixtures new simula- 
tion results are presented. Tables 3 and 4 compile the full simulation data set. 
The statistical uncertainties of this data are estimated to be in the order of 1 
MPa. This translates into a relative error of around 2 % for medium tempera- 
tures, but the relative errors are considerably larger at extremely low and high 
temperatures, where the inversion pressure approaches zero. 

For all systems the results for the JT inversion curves from molecular simulation 
are compared to the four mixture EOS, i.e. DDMIX [18,19], SUPERTRAPP 
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[20], BACKONE-EOS [23] and GERG-2004 [24]. The pure substance results are 
additionally compared to reference EOS using the program package REFPROP 
[50]. All results are presented in pressure-temperature diagrams as well as in 
deviation plots. 



3.1 Pure components 



Figures 1 and 2 present the JT inversion curves of the pure fluids, which were 
grouped to achieve a good visibihty. In these Figures, it can be seen that the 
results from the different methods qualitatively agree. Especially in the low tem- 
perature region of the JT inversion curves, they are often undistinguishable in 
these absolute plots. Significant deviations, however, occur for higher tempera- 
tures. For a more detailed discussion, the deviation plots in Figure 3 are more 
suited, where the baselines represent GERG-2004. It should be pointed out that 
in Figure 3 additionally the pure substance reference EOS [13,14,15,16] results 
arc shown. GERG-2004 agrees with the reference EOS within less than about 3 % 
throughout, except for methane at high temperatures. DDMIX and BAGKONE 
agree roughly equally well to GERG-2004 as well as the simulation data, often 
with the same trends. Generally, these results lie within a band of about 5 %, 
larger deviations are found particularly at high temperatures. BAGKONE yields 
mostly higher results for high temperatures than GERG-2004, whereas DDMIX 
tends to yield lower values. Significantly worse is SUPERTRAPP which yields 
negative deviations of more than 10 % in large parts of the high temperature 
range of methane and nitrogen. 
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3.2 Mixtures 



Figures 4 to 6 present the JT inversion curves for the six equimolar binary and 
ternary mixtures, which are again grouped to achieve good visibihty. As for the 
pure fluids, the quahtative agreement between the different studied models is ob- 
served throughout. The agreement between aU models is almost always excellent 
for low temperatures, but significant deviations are found in the high tempera- 
ture region. Compared to the pure substances the spread between the different 
models is larger. Also here, the results can better be resolved in deviation plots, 
where the baselines were again chosen to represent GERG-2004 data. Figure 7 
shows the binary and Figure 8 the ternary cases. Simulation data and GERG- 
2004 agree also for mixtures usually within 5 %, larger relative deviations are 
found for very high temperatures. BACKONE shows a comparable performance. 
It yields systematically higher values at strongly elevated temperatures and a 
better agreement with GERG-2004 than DDMIX and SUPERTRAPP. These 
two EOS show poorer results throughout, which are always too low by more 
than 10 % at high temperatures. SUPERTRAPP yields the largest deviations in 
all cases. 

GERG-2004 was taken as a reference here as the broadest possible experimental 
data set was taken for its parameterization. This is confirmed by the fact that 
GERG-2004 data lies almost always in between the remaining results. Comparing 
it to the reference EOS for pure fluids, which were individually fltted including 
JT inversion data, an uncertainty of only about 3 % has to be assumed, except 
for methane at high temperatures. 
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4 Conclusion 



In this work, results from molecular modelling and simulation were compared 
to four advanced mixture EOS regarding JT inversion curves of natural gas 
mixtures. With a focus on the four most important components methane, ni- 
trogen, carbon dioxide and ethane, six different equimolar binary and ternary 
mixtures were selected systematically. As a reference and for completeness, the 
four components were also regarded as pure fluids. The comparison shows that 
molecular modelling and simulation is competitive with the most state-of-the-art 
EOS in predicting JT inversion curves. This approach has a similar performance 
as BACKONE, which is based on molecular simulation data itself. DDMIX and 
particularly SUPERTRAPP are less rehable. 

It can be stated that for other mixtures, where no such elaborate EOS are avail- 
able, that molecular modelling and simulation is the method of choice to predict 
JT inversion. 
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List of symbols 



a 


interaction site index 


b 


interaction site index 


Cp 


molar isobaric heat capacity 


fq 


short notation for a trigonometric function 


h 


molar enthalpy 


H 


extensive enthalpy 


ks 


Boltzmann constant 


L 


molecular elongation 


N 


number of molecules 


P 


pressure 


Q 


molecular quadrupole momentum 


Tab 


site-site distance 


rc 


center-center cut-off radius 


T 


temperature 


u 


pair potential 


V 


extensive volume 


X 


mole fraction 


e 


Lennard- Jones energy parameter 




Joule-Thomson coefficient 


i 


binary interaction parameter 


a 


Lennard- Jones size parameter 



Vector properties 

r distance vector 

X mole fraction vector 

U3 orientation vector 
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Subscript 



a interaction site index 

b interaction site index 

i molecule index 

i component index 

j molecule index 

j component index 

Q point quadrupole 

2CLJ two-ccntcr Lcnnard- Jones 



2CLJQ two-center Lennard- Jones plus point quadrupole 
Superscript 
res residual property 
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Tabic 1 

Parameters of pure fluid molecular models, taken from 
[40]. 



Fluid a/A (e/A;B)/K L/A Q/ DA 

methane 3.7281 148.55 

nitrogen 3.3211 34.897 1.0464 1.4397 

carbon dioxide 2.9847 133.22 2.4176 3.7938 

ethane 3.4896 136.99 2.3762 0.8277 



Table 2 

Binary interaction parameters, taken from [41,42,43]. 

Mixture ^ 

methane + carbon dioxide 0.997 

methane + ethane 0.958 

nitrogen + methane 0.974 

nitrogen + carbon dioxide 0.962 

nitrogen + ethane 0.954 

carbon dioxide + ethane 1.041 
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Table 3 

Molecular simulation results for JT inversion of the pure components. Data for methane 
and carbon dioxide was taken from [25], the remainder is from this work. 



T / K 


p / MPa 


T / K 


p / MPa 


T / K 


p 1 MPa 


T / K 


p 1 MPa 


methane 




nitrogen 


carbon 


dioxide 


ethane 




178.26 


15.83 


100 


1.81 


300 


29.05 


250 


2.83 


222.83 


28.43 


125 


12.31 


350 


58.80 


275 


12.50 


267.39 


41.22 


150 


24.68 


400 


76.09 


300 


21.58 


311.96 


47.53 


175 


31.13 


450 


82.32 


325 


29.42 


356.52 


50.83 


200 


35.78 


500 


87.16 


350 


34.94 


401.09 


51.73 


250 


39.54 


550 


90.27 


375 


40.45 


445.65 


51.44 


300 


38.88 


600 


91.31 


425 


48.80 


490.22 


50.58 


350 


36.72 


650 


90.62 


475 


53.93 


534.78 


48.73 


400 


31.55 


700 


87.16 


500 


57.38 


579.35 


45.68 


450 


24.62 


750 


80.94 


525 


58.28 


623.91 


42.81 


500 


17.37 


800 


76.09 


600 


60.14 


668.48 


39.56 


550 


8.48 


850 


69.18 


675 


60.82 


713.04 


33.62 






900 


60.87 


750 


56.97 


757.61 


27.73 






950 


53.96 


825 


54.34 


802.17 


22.68 






1000 


45.66 


900 


48.37 


846.74 


15.27 






1050 


37.35 


975 


42.13 


891.30 


10.68 






1100 


28.36 


1000 


38.86 


935.87 


6.32 






1150 


17.99 






965.58 
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Table 4 

Molecular simulation results for JT inversion of the mixtures, this work. 
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Abstract 

Molecular modelling and simulation as well as four equations of state (EOS) are ap- 
plied to natural gas mixtures regarding Joule-Thomson (JT) inversion. JT inversion 
curves are determined by molecular simulation for six different natural gas mixtures 
consisting of Methane, Nitrogen, Carbon Dioxide and Ethane. These components 
are also regarded as pure fluids, leading to a total of ten studied systems. The results 
are compared to four advanced mixture EOS: DDMIX, SUPERTRAPP, BACKONE 
and the recent GERG-2004 Wide- Range Reference EOS. It is found that molecu- 
lar simulation is competitive with state-of-the-art EOS in predicting JT inversion 
curves. The molecular based approaches (simulation and BACKONE) are superior 
to DDMIX and SUERTRAPP. 

Key words: 



1 Introduction 



Due to the eminent importance of natural gas, knowledge on its thermodynamic 
behaviour and appropriate property models are of great interest. A property often 
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needed for applications of natural gases is the adiabatic (or isenthalpic) Joule- 
Thomson ( JT) coefficient which is defined as the derivative of the temperature 
T with respect to the pressure p at constant enthalpy h and constant composition 

X 



where Cp is the isobaric heat capacity. The JT inversion curve, connecting all state 
points with fi=0, divides the pressure-temperature plane into two regions. In the 
lower region, is positive so that an adiabatic expansion leads to a decrease in 
temperature. In the upper region, /i is negative. It can be shown that the cooling 
effect is maximized if an expansion starts from the inversion pressure. 

The experimental determination of a fluid's JT inversion curve demands precise 
measurements of volumetric and caloric properties over a broad range of tem- 
perature and pressure. Temperature and pressure reach up to a five-fold and 
twelve-fold of their critical values, respectively. Experimental JT inversion curve 
data are therefore scarce, sometimes unreliable [1], and mostly available only for 
pure fluids [2]. An example for mixture data is the work of Charnley et al. [3] 
on e.g. Carbon Dioxide -|- Nitrous Oxide, Carbon Dioxide -|- Ethylene. A good 
overview over experimental data is given in [2]. 

For industrial applications there is a need for a proper representation and also for 
the prediction of JT data, as it is not feasible to measure it for all relevant and 




(1) 



or, using basic thermodynamics relations. 
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often very different blends. Some authors, e.g., Miller [7] or Gunn et al. [8], have 
proposed direct representations of the JT inversion curve, which have hardly any 
predictive power. Significantly more valuable are proper equations of state (EOS) 
that contain much more thermodynamic information and are valid for a broad 
range of state points. Thus, extensive efforts are made to use EOS for predictions 
of the JT inversion curve. Examples for the use of cubic EOS are modified versions 
of Peng-Robinson [9,10], Redlich-Kwong [9,11], Soave-Redlich-Kwong [10], Patel- 
Teja [9] or other cubic EOS [12] with varying parameter functions. Unfortunately, 
cubic EOS are usually unreliable, especially in the high temperature region, where 
they often fail completely. 

Among mixtures, natural gases are the ones that were investigated most exten- 
sively both experimentally and theoretically so that very reliable thermodynamic 
data and models are available. Therefore, natural gases are excellent test systems 
to validate thermodynamic models for mixtures. 

Natural gas from the rig is a mixture of typically seventeen components (con- 
taining Methane, Nitrogen, Carbon Dioxide, Ethane, Propane) [4], but usually 
its main component is Methane [5,6]. As a natural product, it has a great vari- 
ability in composition and, depending on conditions in the formation process, 
considerable quantities of Nitrogen (up to 60 mole %), Carbon Dioxide (up to 
mole 50 %) or Ethane (up to mole 20 %) are encountered [5]. 

For a number of pure natural gas components, reference EOS have been devel- 
oped based on a vast experimental data set considering different thermodynamic 
properties, e.g.. Methane [13], Nitrogen [14], Carbon Dioxide [15] and Ethane 
[16]. Reference EOS have an empirical background, but they are parameterized 
extremely carefully, taking also available experimental JT coefficients into ac- 
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count. Hence, they are regarded in this work as the best available information. 



For mixtures, the National Institute of Standards and Technology (NIST) [17] 
provided two classical phenomenological EOS, i.e. DDMIX [18,19] and SUPER- 
TRAPP [20]. DDMIX is an implementation of the NIST extended corresponding 
states model for mixtures, whereas SUPERTRAPP is based on both a modi- 
fied Peng- Robinson EOS and the NIST extended corresponding states model for 
mixtures. Both were parameterized to experimental pure substance and mixture 
data. Particularly SUPERTRAPP is often used in the literature as a property 
model for designing cooling cycles with mixed coolants, e.g. [21,22]. 

There are also physically based EOS that take the different molecular inter- 
actions, like dispersion or polarity, explicitly into account; an example is the 
BACKONE-EOS [23] . Such EOS can be parameterized for real substances with 
a very small experimental data set, e.g. a few vapour-liquid equilibrium data 
points, as they have a good predictive power. Furthermore, the GERG-2004 
Wide-Range Reference EOS [24] has become available recently, which applies 
the concept of reference EOS to mixtures. Here too, a vast set of experimental 
data, both pure substance and mixture, was used for the development. 

Molecular modelling and simulation offers an interesting alternative approach for 
predicting thermodynamic properties. Instead of describing those macroscopic 
properties directly, the intermolecular interactions are described. Previous work 
from our group [25] has demonstrated the good predictive power of molecular 
models for JT inversion curves for different pure fluids, but also for the mixture 
air which has been modelled as a ternary system containing Nitrogen, Oxygen 
and Argon. The objective of the present paper is to validate the predictive power 
of available molecular models by comparing the results to the four advanced EOS 
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mentioned above using different natural gas mixtures and tfieir most important 
pure components regarding JT inversion. 

2 Molecular model and simulation method 

Most publications on JT inversion curves by molecular simulation are based on 
the spherical Lennard- Jones (LJ) model [26,27,28,29], which is appropriate only 
for simple molecules like Methane and the noble gases. However, this simple po- 
tential model is well suited to further develop molecular simulation techniques for 
determining JT inversion as reliable simulation data is available for comparison. 
E.g., Cohna et al. [30] have chosen two different routes, i.e. via compressibihty 
[31,32] or via thermal expansivity [33], to simulate the JT inversion curve. Work 
on more complex fluids is still scarce, examples are Chacin et al. [34] (Carbon 
Dioxide), Lisal et al. [35] (R32), Kristof et al. [36] (Hydrogen Sulphide), or a 
recent work of our group [25] dealing with fifteen different pure substances (in- 
cluding Methane, Carbon Dioxide, R134a, R143a, R152a, etc.) and the mixture 
air. 

The effective 2CLJQ pair potential was used as molecular model here to describe 
the intermolecular interactions in all cases. This can be done, as only the four 
most common components Methane, Nitrogen, Carbon Dioxide and Ethane were 
considered, being well suited for this modelling approach. The 2CLJQ potential 
is composed out of two identical Lennard- Jones sites a distance L apart (2CLJ) 
and a point quadrupole with momentum Q placed in the geometric centre of the 
molecule oriented along the molecular axis 

U2Ch3q{rij, uji, ujj, L, Q) = 'U2CLj(nj, w^, w^, L) + uq{rij, w^, w^, Q), (3) 
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wherein ii2CLj is the contribution of the four Lennard- Jones interactions 



Herein, is the centre-centre distance vector of two molecules i and j, Tab is one 
of the four Lennard- Jones site-site distances; a counts the two sites of molecule 
i, b counts those of molecule j. The vectors u)i and u)j represent the orientations 
of the two molecules i and j. /q is a trigonometric function depending on these 
molecular orientations, cf . [37] . The Lennard- Jones parameters a and e represent 
size and energy, respectively. 2CLJQ models have the four state independent 
model parameters a, e, L, and Q, which have been adjusted to experimental 
vapour pressure, bubble density and critical temperature in a recent work of our 
group [38]. All pure substance parameters are given in Table 1. The spherical 
non-polar LJ model for Methane is a limiting case of 2CLJQ models, where L — 
and g = 0. 

To perform simulations of mixtures, a molecular mixture model is needed. On 
the basis of pairwise additive pure fluid models, molecular modelling of mixtures 
reduces to modelling interactions between unlike molecules. Here, the modified 
Lorentz-Berthelot combining rule with one adjustable binary interaction param- 
eter ^ was used for each unlike Lennard- Jones interaction 





The quadrupolar contribution uq, is given by [37] 





cta + ctb 



(6) 



ctab — 



2 



^AB — V^A- CB- 



(7) 
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The state independent parameter ^ was adjusted to one experimental vapour 
pressure of each binary mixture in prior work of our group [39,40,41]. Table 2 
reports the six binary parameters of the quaternary natural gas mixture model. 
The unlike quadrupolar interactions are treated in a physically straightforward 
way, following the laws of electrostatics without any binary parameters. It should 
be pointed out, that exclusively experimental VLE data were used in the param- 
eterization of the molecular mixture model, but no caloric data. 

For the calculation of JT inversion curves on the basis of a given model by molecu- 
lar simulation, several methods have been proposed in the literature [31,32,33,42,43]. 
Here, as in prior work [25], an intuitive and straightforward method was used. 
To calculate one JT inversion pressure p for a given temperature T, a series of 
simulations, generally from 5 to 10, were made around the expected result, cov- 
ering typically a rather large pressure range of 20 MPa. In these simulations, the 
enthalpy h and its partial derivative with respect to pressure {dh/dp)T,x were cal- 
culated. Both data sets were fitted simultaneously by a second order polynomial 
vs. pressure at that particular temperature. The inversion pressure corresponds 
simply to the minimum value of enthalpy, i.e. the minimum of that quadratic fit. 

Molecular dynamics simulations were performed in the isobaric-isothermal (NpT) 
ensemble, using Andersen's barostat [44] and isokinetic velocity scahng [45] for 
thermostating. After 6 000 equilibration time steps, the residual enthalpy [46] 

^ ( N N ] 

h""' = M- < EE«2CLjQK-,u;,,a;,-,L,g) > +p < > - fegT, (8) 
( i=ij>i J 

was averaged over 200 000 time steps, where the first term indicates the simula- 
tion average over the intermolecular potential energy and < > is the average 
of the extensive volume. The partial derivative was obtained by a fiuctuation 
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expression [47] 



f dh'^^^\ 1 r 1 1 

\-w),„_ = N- " - ^ >1 + < ^ >} • 

where H'^^^ is the extensive residual enthalpy. The ideal part of the enthalpy is 
not relevant, as it is not pressure dependent. 

Depending on the density of the state point, the membrane mass parameter M 
of Andersen's barostat [44] was chosen from 10~^° to 10~^^ kg/m^. The inter- 
molecular interactions were evaluated explicitly up to a cut-off radius of Tc = 5a 
and standard long range corrections were used for the Lennard- Jones interaction, 
employing angle averaging as proposed by Lustig [48] . Long-range corrections for 
the quadrupolar interaction are not needed since they disappear. A total number 
oi N — 1372 molecules were initially placed in a fee lattice configuration into a 
cubic simulation box, where the density of the system was chosen close to that 
expected from an EOS, if available, otherwise estimates were used. 

3 Results 

JT inversion curves are compared for six different systematically chosen gas mix- 
tures consisting of the four main natural gas components Methane, Nitrogen, 
Carbon Dioxide and Ethane. Firstly, the three Methane containing binary sys- 
tems that can be formed from these four components were studied at equimolar 
composition, where the highest effect of mixing can be expected. Following this, 
secondly, all three ternary systems containing Methane were studied, again at 
equimolar composition, i.e. with a mole fraction Xi — 1/3 throughout. Finally, 
also the four pure fluids were investigated. 
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Simulation results for pure Methane and Carbon Dioxide were taken from prior 
work [25], but for Nitrogen and Ethane as well as for all mixtures new simula- 
tion results are presented. Tables 3 and 4 compile the full simulation data set. 
The statistical uncertainties of this data are estimated to be in the order of 1 
MPa. This translates into a relative error of around 2 % for medium tempera- 
tures, but the relative errors are considerably larger at extremely low and high 
temperatures, where the inversion pressure approaches zero. 

For all systems the results for the JT inversion curves from molecular simulation 
are compared to the four mixture EOS, i.e. DDMIX [18,19], SUPERTRAPP 
[20], BACKONE-EOS [23] and GERG-2004 [24]. The pure substance results are 
additionally compared to reference EOS using the program package REFPROP 
[49]. All results are presented in pressure-temperature diagrams as well as in 
deviation plots. 

3.1 Pure components 

Figures 1 and 2 present the JT inversion curves of the pure fluids, which were 
grouped to achieve a good visibility. In these Figures, it can be seen that the 
results from the different methods qualitatively agree. Especially in the low tem- 
perature region of the JT inversion curves, they are often undistinguishable in 
these absolute plots. Significant deviations, however, occur for higher tempera- 
tures. For a more detailed discussion, the deviation plots in Figure 3 are more 
suited, where the basehnes represent GERG-2004. It should be pointed out that 
in Figure 3 additionally the pure substance reference EOS [13,14,15,16] results 
are shown. GERG-2004 agrees with the reference EOS within less than about 3 % 
throughout, except for Methane at high temperatures. DDMIX and BACKONE 
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agree roughly equally well to GERG-2004 as well as the simulation data, often 
with the same trends. Generally, these results lie within a band of about 5 %, 
larger deviations are found particularly at high temperatures. BACKONE yields 
mostly higher results for high temperatures than GERG-2004, whereas DDMIX 
tends to yield lower values. Significantly worse is SUPERTRAPP which yields 
negative deviations of more than 10 % in large parts of the high temperature 
range of Methane and Nitrogen. 

3.2 Mixtures 

Figures 4 to 6 present the JT inversion curves for the six equimolar binary and 
ternary mixtures, which are again grouped to achieve good visibility. As for the 
pure fluids, the qualitative agreement between the different studied models is ob- 
served throughout. The agreement between all models is almost always excellent 
for low temperatures, but significant deviations are found in the high tempera- 
ture region. Compared to the pure substances the spread between the different 
models is larger. Also here, the results can better be resolved in deviation plots, 
where the baselines were again chosen to represent GERG-2004 data. Figure 7 
shows the binary and Figure 8 the ternary cases. Simulation data and GERG- 
2004 agree also for mixtures usually within 5 %, larger relative deviations are 
found for very high temperatures. BACKONE shows a comparable performance. 
It yields systematically higher values at strongly elevated temperatures and a 
better agreement with GERG-2004 than DDMIX and SUPERTRAPP. These 
two EOS show poorer results throughout, which are always too low by more 
than 10 % at high temperatures. SUPERTRAPP yields the largest deviations in 
all cases. 



10 



GERG-2004 was taken as a reference here as the broadest possible experimental 
data set was taken for its parameterization. This is confirmed by the fact that 
GERG-2004 data lies almost always in between the remaining results. Comparing 
it to the reference EOS for pure fluids, which were individually fitted including 
JT inversion data, an uncertainty of only about 3 % has to be assumed, except 
for methane at high temperatures. 



4 Conclusion 



In this work, results from molecular modelling and simulation were compared to 
four advanced mixture EOS regarding JT inversion curves of natural gas mix- 
tures. With a focus on the four most important components Methane, Nitrogen, 
Carbon Dioxide and Ethane, six different equimolar binary and ternary mix- 
tures were selected systematically. As a reference and for completeness, the four 
components were also regarded as pure fluids. The comparison shows that molec- 
ular modelling and simulation is competitive with the most state-of-the-art EOS 
in predicting JT inversion curves. This approach has a similar performance as 
BACKONE, which is based on molecular simulation data itself. DDMIX and 
particularly SUPERTRAPP are less rehable. 

It can be stated that for other mixtures, where no such elaborate EOS are avail- 
able, that molecular modelling and simulation is the method of choice to predict 
JT inversion. 
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Tabic 1 

Parameters of pure fluid molecular models, taken from 
[38]. 



Fluid a/k {e/k^)/K L/k Q/T)k 

Methane 3.7281 148.55 

Nitrogen 3.3211 34.897 1.0464 1.4397 

Carbon Dioxide 2.9847 133.22 2.4176 3.7938 

Ethane 3.4896 136.99 2.3762 0.8277 



Table 2 

Binary interaction parameters, taken from [39,40,41]. 

Mixture ^ 

Methane + Carbon Dioxide 0.997 

Methane + Ethane 0.958 

Nitrogen + Methane 0.974 

Nitrogen + Carbon Dioxide 0.962 

Nitrogen + Ethane 0.954 

Carbon Dioxide + Ethane 1.041 
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Table 3 

Molecular simulation results for JT inversion of the pure components. Data for Methane 
and Carbon Dioxide was taken from [25] , the remainder is from this work. 



T / K 


p / MPa 


T / K 


p 1 MPa 


T / K 


p 1 MPa 


T / K 


p 1 MPa 


Methane 


Nitrog 


en 


Carbon 


Dioxide 


Ethane 




178.26 


15.83 


100 


1.81 


300 


29.05 


250 


2.83 


222.83 


28.43 


125 


12.31 


350 


58.80 


275 


12.50 


267.39 


41.22 


150 


24.68 


400 


76.09 


300 


21.58 


311.96 


47.53 


175 


31.13 


450 


82.32 


325 


29.42 


356.52 


50.83 


200 


35.78 


500 


87.16 


350 


34.94 


401.09 


51.73 


250 


39.54 


550 


90.27 


375 


40.45 


445.65 


51.44 


300 


38.88 


600 


91.31 


425 


48.80 


490.22 


50.58 


350 


36.72 


650 


90.62 


475 


53.93 


534.78 


48.73 


400 


31.55 


700 


87.16 


500 


57.38 


579.35 


45.68 


450 


24.62 


750 


80.94 


525 


58.28 


623.91 


42.81 


500 


17.37 


800 


76.09 


600 


60.14 


668.48 


39.56 


550 


8.48 


850 


69.18 


675 


60.82 


713.04 


33.62 






900 


60.87 


750 


56.97 


757.61 


27.73 






950 


53.96 


825 


54.34 


802.17 


22.68 






1000 


45.66 


900 


48.37 


846.74 


15.27 






1050 


37.35 


975 


42.13 


891.30 


10.68 






1100 


28.36 


1000 


38.86 


935.87 


6.32 






1150 


17.99 






965.58 


3.17 






1200 


8.30 
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Table 4 

Molecular simulation results for JT inversion of the mixtures, this work. 



r / K 


p 1 MPa 


T / K 


p / MPa 


T / K 


p / MPa 


T / K 


p 1 MPa 


Methane + Carbon Dioxide 


Methane + Ethane 


Nitrogen 


-T ivietnane 


Methane + 


275 


40.30 


275 


29.10 




U.oU 


Carbon Dioxide + 


350 


58.68 


350 


45.14 


zuu 


ol.Zo 


Nitrog 


en 


425 


68.17 


425 


53.33 


/ ( 


A'X AA 


200 


22.87 


500 


70.13 


500 


58.24 


350 


46.36 


250 


39.76 


575 


66.90 


575 


58.26 


425 


42.24 


300 


51.44 


650 


63.22 


625 


57.86 


500 


37.35 


350 


57.41 


725 


55.81 


700 


55.53 


575 


28.46 


425 


59.35 


800 


45.20 


775 


48.55 


650 


17.46 


500 


58.68 


875 


34.86 


850 


40.73 


Methane + 


575 


51.29 


950 


23.07 


925 


33.50 


Carbon Dioxide + 


650 


43.62 


Methane + 




950 


29.86 


Ethane + 


725 


30.74 


Ethane + 




1050 


16.44 


275 


31.62 


800 


21.71 


Nitrogen + 








350 


49.72 


875 


15.03 


200 


17.46 






425 


60.60 






275 


38.53 






500 


65.13 






350 


48.64 






575 


65.15 






425 


52.45 






625 


62.31 






500 


51.80 






700 


59.16 






575 


48.77 






775 


52.53 






650 


42.27 






850 


45.07 






725 


35.46 






925 


35.76 






800 


25.26 






1000 


25.69 
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